E013 Pho4 chimera activity analysis using PHO5 reporter, analysis

require(plotly)
require(tidyverse)
require(ggridges)
require(cowplot)
require(RColorBrewer)
require(grid)
require(ggtext)

old <- theme_set(theme_bw(base_size = 16))

Error in theme_set(theme_bw(base_size = 16)) : 
  could not find function "theme_set"

This is the second part of the analysis. In the first part (see ../input/PHO5-data/20231019-pool-qc-PHO5.Rmd), I did QC and exported the filtered dataset. Here, I will continue working with that dataset and answer our biological questions.

Goal

• Analyze the full chimera set flow results for PHO5pr-mCherry reporter.
• Develop an analysis pipeline to perform QC, correction (if needed) and plotting the results.

Data

Import the background subtracted data

dat0 <- read_tsv("../input/20231023-PHO5-bg-subtracted-data.tsv", col_types = "ccccdddddc")

Filter the data

dat <- filter(dat0, host != "PHO84", flag == "pass", date != "02/10") %>% 
  # based on previous QC, the following sample (both replicates) have high
  # variance - one biological replicate is highly expressed, while the other 
  # two have mNeon, but barely any RFP expression.
  mutate(flag = ifelse(plasmid == "233" & host == "pho2∆", "high.var", flag))

Chimera makeup information

meta <- read_tsv("../input/20230208-chimera-Pho4-makeup.txt", col_types = "ccccc")

Summarize data

Here we would like calculate the ratio of RFP/GFP for each chimera (plasmid) across all replicates, including from different days. Note that the parameter of interest is a ratio, which can be estimated using either “means of ratios” or “ratios of means”. These are just two specific instances of a more general estimator, representing two choices of the weights. The “means of ratios” first calculates the ratios for each replicate within a plasmid, then average them. In this calculation, each replicate is given the weight of 1/n (equal). The “ratios of means” first sum up the GFP and RFP values separately across the replicates for each plasmid, then take the ratio between them. In this estimator, the weight for each replicate is x / sum(x), where x is the denominator in the ratio, i.e., GFP. In other words, this estimator will give more weights to the replicates where the chimera had a higher expression level.

Both estimators are known to be biased. We will ignore that for the moment. In terms of a choice between the two, it seems that there is no reason to give more weights to the experiments with a higher GFP signal. So, the “means of ratios” seems a more natural choice. However, we will calcultae both and dedice later.

A final question is how to calculate the variance of the ratio estimate. According to the survey package manual, an approximate estimator for the variance is

\[ r = \frac{\bar{y}}{\bar{x}}, \text{where}\ \bar{y}=\frac{1}{n}\sum_{i=1}^{n}y_i\ \text{and}\ \bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i\ \\ \hat{V}(r) = (1-\frac{n}{N})(\frac{1}{\bar{x}^2})\frac{s_r^2}{n}\ \text{where}\ s_r^2=\frac{1}{n-1}\sum_{i=1}^{n}(y_i-rx_i)^2 \]

Assuming that N>>n, we can ignore the first term in the variance estimator. The rest can be calculated from the data

datsum <- dat %>%
  filter(!is.na(plasmid)) %>% 
  group_by(plasmid, host) %>% 
  summarize(
     n = n(),
    mG = mean(BL1.H),
    mR = mean(YL2.H),
     A = mean(YL2.H/BL1.H),
     r = mR/mG,
    s2 = 1/(n-1)*sum((YL2.H - r*BL1.H)^2),
    vr = 1/(mG^2)*s2/n,
    se = sqrt(vr),
    .groups = "drop"
  ) %>% 
  select(-s2, -vr)# %>% 
  #pivot_wider(names_from = host, values_from = BL1.H:`nR/G`) %>% 
  #mutate(`pho2∆/PHO2` = `R/G_pho2∆`/`R/G_PHO2`,
  #       `n.pho2∆/PHO2` = `nR/G_pho2∆`/`nR/G_PHO2`)

For each chimera, we would also like to calculate three values:

1. A in pho2∆: this is its base activity without Pho2
2. A in PHO2: this is its full activity with Pho2
3. A_PHO2 / A_pho2∆: this is the Pho2 enhancement of activity

We assign the chimeras into several groups, based on their A_PHO2 and A_PHO2/A_pho2∆

ximera <- datsum %>%
  mutate(host = fct_recode(host, "pho2" = "pho2∆")) %>% 
  pivot_wider(id_cols = plasmid, names_from = host, values_from = A, names_prefix = "A_") %>% 
  mutate(
    s_PHO2 = A_PHO2 / A_PHO2[plasmid == "194"],
    s_pho2 = A_pho2 / A_pho2[plasmid == "194"],
    boost = A_PHO2 / A_pho2,
    group = case_when(
      plasmid %in% c("188", "194") ~ "ref",
      s_PHO2 < 0.2                 ~ "n.f.",
      .default = "chimera"
    ),
    group = fct_relevel(group, "ref", "chimera", "n.f.")
  ) %>% 
  right_join(select(meta, plasmid, set, symbol), by = "plasmid") %>% 
  mutate(symbol = fct_reorder(symbol, s_PHO2, .desc = TRUE)) %>% 
  relocate(c(set, symbol, group), .after = plasmid)

Export the summarized data

write_tsv(ximera, file = "../output/20231125-PHO5pr-chimera-summarized.tsv")

Data selection

Write a function that can select a subset of the chimeras given a set of rules

my_data_select <- function(pattern = NULL, Set = NULL){
  # change region 4 into a consistent format
  tmp <- ximera %>% 
    mutate(
      symbol = as.character(symbol),
      Symbol = ifelse(
        nchar(symbol) == 5,
        paste0(str_sub(symbol, 1, 3), 
               str_sub(symbol, 4, 4), 
               str_sub(symbol, 4, 4), 
               str_sub(symbol, 5, 5)),
        symbol
      )) %>% 
    select(plasmid, Symbol)
  # starting set
  if(length(Set) == 0)
    xim <- filter(tmp, !plasmid %in% refs)
  else
    xim <- filter(tmp, set %in% Set, !plasmid %in% refs)
  # compare to the pattern
  try(if(nchar(pattern) != 6) stop("Pattern must be a string with 6 characters"))
  symbols = xim$Symbol # extrac the symbols for testing
  include = nchar(symbols) > 0 # initialize the inclusion vector
  for(i in 1:6){
    p = substr(pattern, i, i)
    if(p != "X" & p != "x"){ # ignore X and x
      test = toupper(str_sub(symbols, i, i)) == toupper(p)
      include = include & test
    }
  }
  select <- cbind(xim, include)
  return(select$plasmid[select$include])
}
my_data_select_m <- function(patterns = NULL, Set = NULL){
  all_selected = c()
  try(if(length(patterns) == 0) stop("No patterns provided"))
  for(i in patterns)
    all_selected = c(all_selected, my_data_select(pattern = i))
  return(unique(all_selected))
}

Plotting functions

Set up common parameters for thresholding and plotting

# reference Pho4 plasmid ids
refs <- c("188", "194")
# colors
date.colors = c(brewer.pal(name="Dark2", n = 8), brewer.pal(name="Paired", n = 8))
host.colors = c("PHO2" = "gray30", "pho2∆" = "gray70")
point.colors = c("PHO2" = "forestgreen", "pho2∆" = "purple4")
# 

Data prep and transform

my_data_prep <- function(selection){
  # given a selection of chimera ID (plasmid), prepare a data frame for plotting
  # subset data
  tmp <- ximera %>% 
    filter(plasmid %in% c(refs, selection)) %>% 
    select(plasmid, symbol, group) %>% 
    inner_join(dat, by = "plasmid") %>% 
  return(tmp)
}

Plot RFP/GFP ratio and individual components

my_plot_ratio <- function(selection){
  tmp <- my_data_prep(selection)
  p <- tmp %>% 
    select(-c(FSC.H, nGFP, nRFP, flag)) %>% 
    mutate(`R/G` = YL2.H/BL1.H) %>% 
    pivot_longer(cols = c(BL1.H, YL2.H, `R/G`), 
                 names_to = "parameter", values_to = "value") %>% 
    mutate(parameter = factor(parameter, levels = c("R/G", "YL2.H", "BL1.H"),
                              labels = c("RFP/GFP", "PHO5pRFP", "Pho4-GFP"))) %>% 
    ggplot(aes(x = symbol, y = value, group = host)) + 
    stat_summary(aes(group = host), fun.data = "mean_cl_boot", geom = "errorbar",
                 position = position_dodge(0.5), width = 0.3) +
    geom_bar(aes(fill = host), width = 0.5, alpha = 0.8,
             stat = "summary", fun = "mean", position = position_dodge(0.5)) +
    geom_point(data = function(x) subset(x, !symbol %in% c("CCCCC", "SSSSS")),
               aes(group = host, color = host), size = 1, shape = 3, alpha = 0.9,
               position = position_jitterdodge(dodge.width = 0.5, jitter.width = 0.1)) +
    scale_color_manual(values = point.colors) +
    scale_fill_manual(values = host.colors) +
    facet_grid(parameter~group, scales = "free", space = "free_x") +
    theme_bw(base_size = 18) + background_grid(minor = "none") + 
    xlab("Pho4 chimera") +
    theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"),
          legend.position = "top",
          axis.title = element_blank())
  return(p) 
}

Plot chimera relative activity to ScPho4 w/ Pho2 and boost factor

my_plot_rel_act <- function(selection){
  # given a selection of chimera IDs, plot their functionality w/PHO2
  # relative to ScPho4, and their boost
  p <- filter(ximera, plasmid %in% c(refs, selection)) %>%
    ggplot(aes(x = symbol, y = s_PHO2)) +
    geom_col(width = 0.3, color = "black", fill = "gray80") + 
    geom_hline(yintercept = 1, linetype = 2, color = "gray30") +
    facet_grid(.~group, scales = "free_x", space = "free_x") + 
    scale_y_continuous(labels = scales::percent) +
    xlab("Pho4 chimera") + ylab("A<sub>PHO2</sub>, chimera/ScPho4") +
    ggtitle("Chimera function in <em>S. cerevisiae</em>") +
    theme_bw(base_size = 18) + background_grid(minor = "none") +
    theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"),
          #strip.text = element_blank(), 
          plot.title = element_markdown(hjust = 0.5),
          axis.title.y = element_markdown(),
          legend.position = "top")
  return(p)
}

my_plot_boost <- function(selection){
  # given a selection of chimera IDs, plot their functionality w/PHO2
  # relative to ScPho4, and their boost
  # dat
  tmp <- filter(ximera, plasmid %in% c(refs, selection)) %>% 
    mutate(perc_pho2 = A_pho2/A_PHO2) %>% 
    pivot_longer(cols = c(s_PHO2, boost, perc_pho2), 
                 names_to = "parameter", values_to = "ratio")
  # labeller
  par.explain <- c(
    s_PHO2 = "Rel. A<sub>PHO2</sub>",
    boost = "Boost",
    perc_pho2 = "%A<sub>pho2∆</sub>"
  )
  p <- ggplot(tmp, aes(x = symbol, y = ratio)) +
    geom_col(width = 0.3, color = "black", fill = "gray80") +
    geom_hline(yintercept = 1, linetype = 2, color = "gray30") +
    facet_grid(parameter~group, scales = "free", space = "free_x",
              labeller = labeller(parameter = par.explain)) +
    theme_bw(base_size = 18) +
    background_grid(minor = "none") +
    theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"),
          axis.title.x = element_blank(),
          strip.text.y = element_markdown(size = rel(0.9))
          #plot.title = element_text(hjust = 0.5)
    )
  return(p)
}

High variance samples

Summarize the background subtracted data by calculating the means and cv for each strain.

cv <- dat %>% 
  select(-nGFP, -nRFP) %>%
  pivot_longer(FSC.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  group_by(date, plasmid, host, parameter) %>% 
  summarize(
    n = n(),
    mean = mean(intensity),
    cv = sd(intensity)/mean(intensity),
    .groups = "drop"
  ) %>% 
  arrange(desc(cv))

Number of replicates left for each sample

expt <- dat %>% 
  filter(host %in% c("PHO2", "pho2∆"), !plasmid %in% c("188", "194")) %>% 
  group_by(date, plasmid, host) %>% 
  summarize(n = n(), .groups = "drop")

expt %>% 
  ggplot(aes(x = plasmid, y = n)) +
  geom_col(aes(fill = host)) + 
  facet_grid(date ~ .) +
  scale_fill_manual(values = c("PHO2" = "gray30", "pho2∆" = "gray70")) +
  theme_minimal() + background_grid(major = "none") + panel_border(size = 0.5) +
  scale_y_continuous(name = "Replicates", breaks = c(6)) + xlab(NULL) +
  theme(axis.text.x = element_text(angle = 90),
        strip.text.y = element_text(angle = 0),
        legend.position = "top")

Use the control strain (pH194 with PHO2) to identify and correct for systematic biases

control <- filter(dat, plasmid == "194", host == "PHO2") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  droplevels()

Model for mNeon

gfp.model.0 <- lm(BL1.H ~ log10(events) + date + row*col, data = control)
step(gfp.model.0)

Start:  AIC=1507.62
BL1.H ~ log10(events) + date + row * col

                Df Sum of Sq     RSS    AIC
- row:col        6     96964 3731743 1499.4
- log10(events)  1     12751 3647531 1506.1
<none>                       3634779 1507.6
- date          11   3415318 7050097 1581.0

Step:  AIC=1499.41
BL1.H ~ log10(events) + date + row + col

                Df Sum of Sq     RSS    AIC
- row            3     92309 3824052 1496.9
- log10(events)  1      9919 3741662 1497.8
<none>                       3731743 1499.4
- col            2   1030496 4762240 1530.5
- date          11   3409784 7141528 1570.9

Step:  AIC=1496.93
BL1.H ~ log10(events) + date + col

                Df Sum of Sq     RSS    AIC
- log10(events)  1     16207 3840259 1495.5
<none>                       3824052 1496.9
- col            2   1024530 4848582 1527.1
- date          11   3423810 7247862 1567.0

Step:  AIC=1495.54
BL1.H ~ date + col

       Df Sum of Sq     RSS    AIC
<none>              3840259 1495.5
- col   2   1068825 4909084 1526.9
- date 11   3422018 7262277 1565.3

Call:
lm(formula = BL1.H ~ date + col, data = control)

Coefficients:
(Intercept)    date02/09    date02/11    date02/16    date02/18    date02/19    date02/20    date02/21    date02/22    date02/23    date03/30    date03/31         col5  
     1827.3       -209.7       -259.9       -351.5       -242.3       -319.0       -499.0       -343.8       -233.1       -439.5       -533.6       -569.1       -107.4  
       col9  
     -211.0  

gfp.model.1 <- lm(BL1.H ~ date + col, data = control)

Model for PHO5pr::RFP

rfp.model.0 <- lm(YL2.H ~ log10(events) + date + row*col, data = control)
step(rfp.model.0)

Start:  AIC=2229.67
YL2.H ~ log10(events) + date + row * col

                Df Sum of Sq        RSS    AIC
- row:col        6  13201086  560400308 2221.1
<none>                        547199222 2229.7
- log10(events)  1  13838563  561037785 2231.3
- date          11 852447575 1399646797 2342.9

Step:  AIC=2221.11
YL2.H ~ log10(events) + date + row + col

                Df Sum of Sq        RSS    AIC
<none>                        560400308 2221.1
- log10(events)  1  12713396  573113704 2222.3
- col            2  65908721  626309029 2233.1
- row            3 154001100  714401408 2250.1
- date          11 851423159 1411823467 2332.2

Call:
lm(formula = YL2.H ~ log10(events) + date + row + col, data = control)

Coefficients:
  (Intercept)  log10(events)      date02/09      date02/11      date02/16      date02/18      date02/19      date02/20      date02/21      date02/22      date02/23  
      19089.8         3496.2        -6513.6        -4709.4        -6221.7        -3241.9        -6452.8        -4558.0        -4776.0        -3995.3        -5054.2  
    date03/30      date03/31           rowC           rowE           rowG           col5           col9  
      -8507.9       -10061.8        -1882.5        -2429.8        -2641.7         -954.5        -1663.5  

rfp.model.1 <- lm(YL2.H ~ log10(events) + date + row + col, data = control)

there are more systematic shifts in the RFP, significant for row, col, date and also # of events however, I won’t be removing these effects yet, because I’ve found that RFP/GFP ratios are pretty consistent across days. In other words, the variation in GFP and RFP may be cancelled out.

Check for each plasmid how consistent are the measurements between days

tmp <- dat %>% 
  # remove one sample with only one valid day of experiment
  filter(!(plasmid == "218" & host == "PHO2"), !plasmid %in% c("188", "194", NA)) %>% 
  nest(data = c(date, BL1.H, YL2.H), .by = c(plasmid, host))

day.var.gfp <- tmp %>% 
  mutate(model = map(data, function(df) lm(BL1.H ~ date, data = df)),
         tidied = map(model, broom::tidy)) %>% 
  unnest(tidied) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>% 
  select(-data, -model) %>% 
  filter(p.adj < 0.10) %>% 
  arrange(plasmid, host)

day.var.rfp <- tmp %>% 
  mutate(model = map(data, function(df) lm(YL2.H ~ date, data = df)),
         tidied = map(model, broom::tidy)) %>% 
  unnest(tidied) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>% 
  select(-data, -model) %>% 
  filter(p.adj < 0.10) %>% 
  arrange(plasmid, host)

# extract ximera names
refs <- c("188","194")
# make a test set
day.var.gfp.list <- unique(day.var.gfp$plasmid)
day.var.rfp.list <- unique(day.var.rfp$plasmid)

High day-to-day GFP variance: 212, 222, 229, 231, 241, 251, 252, 277, 301, 326, 328, 329, 331, 334 High day-to-day RFP variance: 212, 216, 239, 241

Plotting components for chimeras with high day-to-day variance in Pho4-mNeon

p <- my_plot_ratio(c(refs,day.var.gfp.list)) + 
    geom_point(data = function(x) subset(x, !symbol %in% c("CCCCC", "SSSSS")),
               aes(group = host, color = date), size = 1, shape = 3, alpha = 0.9,
               position = position_jitterdodge(dodge.width = 0.5, jitter.width = 0.1)) +
    scale_color_manual(values = date.colors, guide = "none")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

p

Watch out for CSCscC, SCCsS, SCCsS

Plotting components for chimeras with high day-to-day variance in PHO5pr-mCherry

p <- my_plot_ratio(c(refs,day.var.rfp.list)) + 
    geom_point(data = function(x) subset(x, !symbol %in% c("CCCCC", "SSSSS")),
               aes(group = host, color = date), size = 1, shape = 3, alpha = 0.9,
               position = position_dodge(width = 0.5)) +
    scale_color_manual(values = date.colors, guide = "none")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

p

most of the day-to-day variance are canceled out after RFP/GFP normalization

All chimera, scatter plot

Design the plot

my_scatter_plot <- function(pattern){
  selection = my_data_select(pattern = pattern)
  scatter.colors = c("ScPho4" = "forestgreen", "CgPho4" = "blue3", "cyan2", "other" = "gray20")
  names(scatter.colors)[3] = pattern
  p <- ximera %>% 
    mutate(A_PHO2 = signif(A_PHO2, digits = 2),
           A_pho2 = signif(A_pho2, digits = 2),
           group = case_when(
             symbol == "CCCCC" ~ "CgPho4",
             symbol == "SSSSS" ~ "ScPho4",
             plasmid %in% selection ~ pattern,
             .default = "other"
           ),
           group = fct_relevel(group, names(scatter.colors))) %>% 
    ggplot(aes(x = A_PHO2, y = A_pho2, label = symbol)) + 
    geom_point(aes(color = group), size = 2.5) + 
    scale_color_manual("Pho4 type", values = scatter.colors) +
    geom_abline(slope = 1) +
    theme_gray(base_size = 14)
  return(p)
}

my_scatter_plot_fix <- function(){
  # this function is the same as above, but is used to plot region 4
  # effects alone, and doesn't take any input
  s1 = my_data_select(pattern = "XXXCCX")
  s2 = my_data_select(pattern = "XXXSSX")
  scatter.colors = c("ScPho4" = "forestgreen", "CgPho4" = "blue3", 
                     "P2ID:Cg" = "cyan2", "P2ID:Sc" = "palegreen",
                     "P2ID:mixed" = "gray20")
  p <- ximera %>% 
    mutate(A_PHO2 = signif(A_PHO2, digits = 2),
           A_pho2 = signif(A_pho2, digits = 2),
           group = case_when(
             symbol == "CCCCC" ~ "CgPho4",
             symbol == "SSSSS" ~ "ScPho4",
             plasmid %in% s1 ~ "P2ID:Cg",
             plasmid %in% s2 ~ "P2ID:Sc",
             .default = "P2ID:mixed"
           ),
           group = fct_relevel(group, names(scatter.colors))) %>% 
    ggplot(aes(x = A_PHO2, y = A_pho2, label = symbol)) + 
    geom_abline(slope = 1) +
    geom_point(aes(color = group), size = 2.5) + 
    scale_color_manual(NULL, values = scatter.colors) +
    labs(x = bquote(A[PHO2]), y = bquote(A[pho2])) +
    theme_minimal(base_size = 16) +
    theme(legend.text = element_text(size = rel(0.75)))
  return(p)
}

p <- my_scatter_plot_fix()
ggsave(filename = "../img/20231220-all-chimera-scatter-color-by-P2ID.png",
       plot = p, width = 6, height = 4, dpi = 300)
ggplotly(p + labs(x = "A<sub>PHO2</sub>", y = "A<sub>pho2</sub>") +
           theme_gray(base_size = 16) +
           theme(legend.text = element_markdown()), 
         tooltip = c("label", "x", "y"))

ANOVA

split <- c(1,1,1,1,1); names(split) <- paste0("P", 1:5)
tmp <- ximera %>% 
  filter(set == "M", group != "n.f.") %>% 
  separate_wider_position(symbol, split) %>% 
  mutate(across(P1:P5, ~factor(.x, levels = c("S", "C"))))
lm <- lm(A_pho2 ~ (P1+P2+P3+P4+P5), data = tmp)
summary(lm)

Call:
lm(formula = A_pho2 ~ (P1 + P2 + P3 + P4 + P5), data = tmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0964 -0.7844 -0.1948  0.5021  5.2278 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.4592     0.9002  -0.510  0.61695    
P1C           2.3355     0.7457   3.132  0.00643 ** 
P2C           1.2950     0.7286   1.778  0.09449 .  
P3C           2.8751     0.7457   3.856  0.00140 ** 
P4C           5.6326     0.9414   5.983 1.91e-05 ***
P5C          -0.1689     0.8543  -0.198  0.84573    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.709 on 16 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.8443,    Adjusted R-squared:  0.7956 
F-statistic: 17.35 on 5 and 16 DF,  p-value: 5.85e-06

The main effects were calculated by averaging over all chimeras with CgPho4 region at the respective position. I’d like to break them down by backgrounds. For example, for region 3, I’d like to see the pairwise comparisons between CCCSS and CCSSS, where only region 3 differs. The steps are

1. select the region to be compared. split the symbol into two parts - the genotype of the focal region and the rest
2. group by the second part (rest) and calculate the differential

my_calc_region_effect <- function(region, variable){
  # this function takes the name of a variable of interest
  # x specifies the foreground region, which will be examined for its effect on
  # the variable of interest.
  # it then transforms the ximera data frame to preserve only the variable of
  # interest, pivots it wider after grouping by the background composition.
  
  # prepare the data by mutating the symbol column into fg and bg
  valid.var <- c("A_PHO2", "A_pho2", "s_PHO2", "spho2", "boost")
  if(!variable %in% valid.var)
    stop(paste0("Please specify one of the valid variable names:", 
                paste(valid.var, collapse = ", ")))
  tmp <- ximera %>% 
    filter(set == "M") %>% 
    select(plasmid, symbol, var = {{ variable }}) %>% 
    mutate(fg = str_sub(symbol, region, region) %>% toupper(),
           bg = symbol %>% toupper())
  # replace the foreground region with X for grouping
  str_sub(tmp$bg, region, region) <- "X"
  # reorganize the tibble for easier handling, optional
  tmp <- relocate(tmp, fg, bg, .before = symbol) %>% select(-symbol)
  # pivot the data into a wide format such that for each background, there
  # are two values for the variable of interest, one from the chimera with 
  # CgPho4's version in the foreground and another with ScPho4's version
  tmp <- tmp %>% 
    select(plasmid, fg, bg, var) %>% 
    pivot_wider(id_cols = bg, names_from = "fg", 
                values_from = c(plasmid, var)) %>% 
    unite(plasmid, starts_with("plasmid")) %>%
    mutate(label = paste(bg, plasmid, sep = "\n"))
  return(tmp)
}

my_plot_region_effect_onevar <- function(region, variable){
  # this function uses `my_calc_region_effect` output as the data
  # and makes a xy scatter plot, where x shows the value of the variable of 
  # interest with CgPho4 in the focal region, and y for the ScPho4 version
  tmp <- my_calc_region_effect(region, variable)
  p <- ggplot(tmp, aes(x = var_C, y = var_S, label = label)) +
    geom_point(size = 2.5) + 
    geom_abline(slope = 1) +
    xlab(paste0("Region ", region, " from CgPho4")) +
    ylab(paste0("Region ", region, " from ScPho4")) +
    xlim(0, NA) + ylim(0, NA) +
    ggtitle(paste0("Effect on ", variable)) +
    theme_gray(base_size = 16) +
    theme(plot.title = element_text(hjust = 0.5))
  return(p)
}

x = 5
p1 <- my_plot_region_effect_onevar(x, "A_PHO2")
p2 <- my_plot_region_effect_onevar(x, "A_pho2")
subplot(p1, p2, margin = 0.05) %>% 
  layout(title = paste("Region", x, "swap effect on A_PHO2 and A_pho2", sep = " "),
         xaxis = list(title = paste0("Region ", x, " from CgPho4")),
         yaxis = list(title = paste0("Region ", x, " from ScPho4")) )

Here, I’d like to take what I build above and create a new tibble, in which each row is a different background (makeup of the chimera except for the focal region). The value columns are:

1. dA_PHO2 = A_PHO2_Cg - A_PHO2_Sc
2. dA_pho2 = A_pho2_Cg - A_pho2_Sc
3. A_PHO2_Sc = A_PHO2_Sc

The goal is to plot dA_PHO2 and dA_pho2 side-by-side for each background.

my_comp_region_effect <- function(region){
  # this function uses my_calc_region_effect to get the value for the variable of interest
  # with either Cg or Sc version in the focal region, separately for each background composition
  # it does so for two variables, A_PHO2 and A_pho2, then calculate dA_PHO2, dA_pho2, and
  # combine them
  PHO2 = my_calc_region_effect(region, "A_PHO2") %>% 
    mutate(dA_PHO2 = var_C - var_S,
           # mean A_PHO2
           M_PHO2 = (var_S + var_C)/2,
           NF = ifelse(M_PHO2 <=3.5, TRUE, FALSE)) %>% 
    select(-var_S, -var_C)
  
  pho2 = my_calc_region_effect(region, "A_pho2") %>% 
    mutate(dA_pho2 = var_C - var_S, 
           M_pho2 = (var_S + var_C)/2) %>% 
    select(-var_S, -var_C)
  
  dat <- full_join(PHO2, pho2, by = c("bg", "plasmid", "label")) %>% 
    select(bg, plasmid, dA_PHO2, dA_pho2, M_PHO2, M_pho2, NF)
  
  return(dat)
}

my_plot_region_effect_twovar_line <- function(region, highlight = "none"){
  # this function uses my_comp_region_effect to generate the data
  # and plot the difference in A_PHO2 and A_pho2 between the CgPho4 vs ScPho4
  # in the focal region
  dat <- my_comp_region_effect(region) %>% 
    pivot_longer(cols = c(dA_PHO2, dA_pho2), 
                 names_to = "host", values_to = "diff") %>% 
    mutate(host = fct_recode(host, `PHO2` = "dA_PHO2", `pho2∆` = "dA_pho2"),
           host = fct_relevel(host, "PHO2"))
  if(highlight != "none" & highlight != region){
    hl = as.numeric(highlight)
    dat <- mutate(dat, grp = str_sub(bg, hl, hl) %>% toupper(),
                  grp = fct_recode(grp, CgPho4 = "C", ScPho4 = "S"))
  }else{
    dat <- mutate(dat, grp = ifelse(NF, "n.f.", "others"))
  }
  # specify legend title
  legend.title = ""
  if(highlight != "none" & highlight != region){
    hl = as.numeric(highlight)
    dat <- mutate(dat, grp = str_sub(bg, hl, hl) %>% toupper(),
                  grp = fct_recode(grp, CgPho4 = "C", ScPho4 = "S"))
    legend.title = paste("Region", highlight, sep = " ")
  }else{
    dat <- mutate(dat, grp = ifelse(NF, "no", "yes"))
    legend.title = "Functional"
  }
  # specify arrow annotation
  arrow.x = 0.7
  arrow.y = (max(dat$diff) - min(dat$diff)) / 5 
  p <- dat %>% 
    ggplot(aes(x = host, y = diff, label = bg)) +
    geom_point(aes(color = grp), size = 2, alpha = 0.8,
               position = position_jitter(0.05)) + 
    geom_line(aes(group = bg), linewidth = 0.2, alpha = 0.8) +
    geom_segment(aes(x = arrow.x, xend = arrow.x, y = -arrow.y, yend = arrow.y),
                 arrow = arrow(length = unit(0.03, "npc"), ends = "both"), 
                 color = "gray60", lwd = 1, alpha = 0.5) +
    geom_segment(aes(x = arrow.x - 0.05, xend = arrow.x + 0.05, y = 0, yend = 0),
                 lwd = 2, color = "gray60") +
    annotate("text", x = arrow.x - 0.1, y = 5, label = "CgPho4++", 
             angle = '90', color = "gray30") +
    annotate("text", x = arrow.x + 0.1, y = -5, label = "ScPho4++", 
             angle = '270', color = "gray30") +
    scale_color_manual(legend.title, values = c("orange", "gray20")) +
    ylab("Region swap effect (Cg-Sc)") +
    theme_bw(base_size = 18) + 
    theme(
      axis.title.x = element_blank(),
      axis.title.y = element_text(size = rel(0.9)),
      legend.text = element_text(size = rel(0.8)),
      legend.title = element_text(size = rel(0.9)),
    )
  return(p)
}

my_plot_region_effect_twovar_line("4", "5")# %>% ggplotly()
ggsave("../img/20231221-region-swap-effect-4-on-5.png", width = 6, height = 4, dpi = 150)

my_plot_region_effect_twovar_line("5", "4")# %>% ggplotly()
ggsave("../img/20231224-region-swap-effect-5-on-4.png", width = 6, height = 4, dpi = 200)

my_plot_region_effect_twovar_line_par <- function(regions, highlight = "none"){
  # this function uses my_comp_region_effect to generate the data
  # and plot the difference in A_PHO2 and A_pho2 between the CgPho4 vs ScPho4
  # in the focal region
  dat <- map_dfr(regions, \(region) my_comp_region_effect(region), .id = "region") %>% 
    pivot_longer(cols = c(dA_PHO2, dA_pho2), 
                 names_to = "host", values_to = "diff") %>% 
    mutate(host = fct_recode(host, `PHO2` = "dA_PHO2", `pho2∆` = "dA_pho2"),
           host = fct_relevel(host, "PHO2"))
  # specify legend title
  legend.title = ""
  if(highlight != "none" & !highlight %in% regions){
    hl = as.numeric(highlight)
    dat <- mutate(dat, grp = str_sub(bg, hl, hl) %>% toupper(),
                  grp = fct_recode(grp, CgPho4 = "C", ScPho4 = "S"))
    legend.title = paste("Region", highlight, "from", sep = " ")
  }else{
    dat <- mutate(dat, grp = ifelse(NF, "no", "yes"))
    legend.title = "Functional"
  }
  # specify arrow annotation
  arrow.x = 0.7
  arrow.y = (max(dat$diff) - min(dat$diff)) / 5 
  p <- dat %>% 
    ggplot(aes(x = host, y = diff, label = bg)) +
    geom_point(aes(color = grp), size = 2, alpha = 0.8,
               position = position_jitter(0.1)) + 
    geom_line(aes(group = bg), linewidth = 0.2, alpha = 0.8) +
    facet_grid(grp ~ region, labeller = labeller(
      grp = c(CgPho4 = "P2ID:CgPho4", ScPho4 = "P2ID:ScPho4"),
      region = label_both
    )) +
    scale_color_manual("P2ID:", values = c("orange", "gray20")) +
    ylab("Region swap effect (Cg-Sc)") +
    theme_bw(base_size = 18) + 
    theme(
      axis.title.x = element_blank(),
      axis.title.y = element_text(size = rel(0.9)),
      legend.text = element_text(size = rel(0.8)),
      legend.title = element_text(size = rel(0.9)),
    )
  return(p)
}
my_plot_region_effect_twovar_line_par(c(1,2,3), "4")
ggsave("../img/20231224-region-swap-effect-1to3-on-4.png")

Saving 7.29 x 4.51 in image

my_plot_region_effect_twovar_side <- function(region){
  # this function uses my_comp_region_effect to generate the data
  # and plot the difference in A_PHO2 and A_pho2 between the CgPho4 vs ScPho4
  # in the focal region
  dat <- my_comp_region_effect(region) %>% 
    pivot_longer(cols = c(dA_PHO2, dA_pho2), 
                 names_to = "host", values_to = "diff") %>% 
    mutate(host = fct_recode(host, `PHO2` = "dA_PHO2", `pho2∆` = "dA_pho2"),
           host = fct_relevel(host, "PHO2"))
  p <- dat %>% 
    ggplot(aes(x = bg, y = diff, group = host)) +
    geom_col(aes(fill = host), position = position_dodge(0.9)) +
    scale_fill_manual(values = host.colors) +
    ylab("Region swap diff (Cg vs Sc)") +
    theme_cowplot(font_size = 20) + 
    panel_border(color = "gray30") +
    background_grid(major = "y", minor = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, family = "courier"),
          axis.title.x = element_blank(),
          legend.position = "top")
  return(p)
}

triangle heatmap

First, write a function to generate the data for plotting. If we are going to use ggplot, we need a tibble to store the data, something in the following form

plasmid	symbol	RegionA	RegionB	A_PHO2	A_pho2	s_PHO2	boost	perc_pho2
209	CCSCC	3	3	8.25	7.82	0.468	1.06	0.94

If we are ok with using non ggplot - heatmaps are not ggplot’s strength anyways - we can just build a matrix.

Note that this way of summarizing the data has many limitaitons: 1) it requires specifying the reference, either CCCCC or SSSSS. Everything is measured against that; 2) it only shows pairwise (two region) interactions. This turns out to be fine with five regions, since every chimera can be expressed as either a 0, 1 or 2 region swap from one of the two reference genotypes. With 6 or more regions, higher level (3 or more region) interactions cannot be visualized this way. Because of this, we will focus on just the main set for this analysis.

To build the matrix, we need to first identify the chimeras that belong to the set. For that, we will use the “main” set, with the five region split, for the moment at least. The function will first determine which reference to use. If we use SSSSS as the reference, for example, we will assign 0 to the reference. All other chimeras with 1 or 2 regions from Cg will be used to fill an upper triangular matrix, using one of the values of interest, e.g., A_PHO2.

my_upper_triangular_mat <- function(alt = "C", var = "A_PHO2"){
  # given the alternative allele (C/S) and a variable of interest, e.g., A_PHO2,
  # output an upper triangular matrix containing the values from the variable 
  # of interest, with the row and col numbers based on the first and second
  # positions containing the alternative allele. If all positions contain the 
  # reference allele, the value is subtracted from all values in the matrix
  # when just one position is the alternative allele, the value in the diagonal
  # is set. when there are more than 2 regions containing the alternative allele
  # skip.
  out_mat <- matrix(NA, nrow = 5, ncol = 5)
  ref_val <- NA
  dat <- filter(ximera, set == "M") %>% 
    mutate(S = as.character(symbol) %>% toupper())
  for(i in seq(1, nrow(dat))){
    symbol = dat[i, "S"]
    # determine which positions contain the alternative allele
    p = str_locate_all(symbol, alt)[[1]][,"start"]
    l = length(p)   # how many positions contain the alt allele
    v = dat[[var]][i] # retrieve the value of the variable
    if(l == 0)
      ref_val = v
    else if(l == 1)
      out_mat[p, p] = v
    else if(l == 2)
      out_mat[p[1], p[2]] = v
  }
  out_mat = out_mat - ref_val
  return(out_mat)
}

my_combined_triangular_mat <- function(alt = "C"){
  # given the alternative allele (C/S), output a matrix containing the values
  # for both with and without Pho2, arranged in two complementary triagular
  # matrices, with the row and col numbers based on the first and second
  # positions containing the alternative allele. If all positions contain the 
  # reference allele, the value is subtracted from all values in the matrix
  # when just one position is the alternative allele, the value in the diagonal
  # is set. when there are more than 2 regions containing the alternative allele
  # skip.
  out_mat <- matrix(NA, nrow = 6, ncol = 6)
  upper <- cbind(NA, my_upper_triangular_mat(alt, var = "A_PHO2")) %>% 
    rbind(., NA)
  lower <- rbind(NA, t(my_upper_triangular_mat(alt, var = "A_pho2"))) %>% 
    cbind(., NA)
  out_mat = ifelse(is.na(upper), lower, upper)
  return(out_mat)
}

my_plot_triangle_heatmap <- function(alt, var){
  # this function takes the output of the function above and makes a heatmap
  # using pheatmap function, then rotates it using grid graphics
  # thanks to https://bookdown.org/rdpeng/RProgDA/the-grid-package.html#grid-graphics-coordinate-systems
  # adding title based on https://davetang.github.io/muse/pheatmap.html
  
  # construct title of plot
  ref = ifelse(alt == "C", "ScPho4", "CgPho4")
  bg = ifelse(var == "A_PHO2", "with PHO2", "w/o pho2")
  my_title <- paste("Epistasis between regions on", ref, "background", bg)
  test <- my_upper_triangular_mat(alt = alt, var = var)
  paletteLength = 50
  myColors <- colorRampPalette(c("steelblue4", "gray90", "red"))(paletteLength)
  rng <- max(abs(test), na.rm = TRUE)
  myBreaks <- c(seq(-rng, 0, length.out=ceiling(paletteLength/2) + 1), 
                seq(rng/paletteLength, rng,
                    length.out=floor(paletteLength/2)))
  p <- pheatmap::pheatmap(test, color = myColors, breaks = myBreaks,
                          border_color = NA, na_col = NA, silent = TRUE,
                          cluster_cols = FALSE, cluster_rows = FALSE)
  vp <- viewport(x = 0.5, y = 0.25,
                 width = unit(4.5, "in"), height = unit(4.5, "in"), angle = 47) 
  grid.newpage()
  pushViewport(vp)
  grid.draw(p$gtable)
  popViewport()
  grid.text(label = my_title, x = 0.5, y = 0.95, gp = gpar(fontsize = 16, fontface = "bold"))
  return(p)
}

my_plot_combined_triangle_heatmap <- function(alt){
  # this function takes the output of the function my_combined_triangular_mat()
  # using pheatmap function, then rotates it using grid graphics
  # thanks to https://bookdown.org/rdpeng/RProgDA/the-grid-package.html#grid-graphics-coordinate-systems
  # adding title based on https://davetang.github.io/muse/pheatmap.html
  
  # construct title of plot
  ref = ifelse(alt == "C", "ScPho4", "CgPho4")
  my_title <- paste("Epistasis between regions on", ref, "background")
  test <- my_combined_triangular_mat(alt = alt)
  paletteLength = 50
  myColors <- colorRampPalette(c("steelblue4", "gray90", "red"))(paletteLength)
  rng <- max(abs(test), na.rm = TRUE)
  myBreaks <- c(seq(-rng, 0, length.out=ceiling(paletteLength/2) + 1), 
                seq(rng/paletteLength, rng,
                    length.out=floor(paletteLength/2)))
  p <- pheatmap::pheatmap(test, color = myColors, breaks = myBreaks,
                          border_color = NA, na_col = NA, silent = TRUE,
                          cluster_cols = FALSE, cluster_rows = FALSE)
  vp <- viewport(x = 0.5, y = 0.45,
                 width = unit(3, "in"), height = unit(2.8, "in"), angle = 47) 
  grid.newpage()
  pushViewport(vp)
  grid.draw(p$gtable)
  popViewport()
  grid.text(label = my_title, x = 0.5, y = 0.95, 
            gp = gpar(fontsize = 16, fontface = "bold"))
  grid.text(label = "With Pho2", x = 0.1, y = 0.65, just = c("left", "top"),
            gp = gpar(fontsize = 14, fontface = "bold"))
  grid.text(label = "Without pho2", x = 0.1, y = 0.25, just = c("left", "top"), 
            gp = gpar(fontsize = 14, fontface = "bold"))
  return(p)
}

p1 <- my_plot_combined_triangle_heatmap("C")

p2 <- my_plot_combined_triangle_heatmap("S")

p1 <- my_plot_triangle_heatmap("C", "A_PHO2")
p2 <- my_plot_triangle_heatmap("C", "A_pho2")
p3 <- my_plot_triangle_heatmap("S", "A_PHO2")
p4 <- my_plot_triangle_heatmap("S", "A_pho2")